##################################################################
##
## Conflict in Lake Chad region
##
## B. Blankespoor 2020-11-05 (first version 2020-04-01)
##
##################################################################


##
## IMPORTANT PATH to Rtools to run Cpp
##
Sys.setenv(PATH = paste("C://WBG//Utility//Rtools//bin", Sys.getenv("PATH"), sep=";"))
Sys.setenv(BINPREF = "C://WBG//Utility//Rtools//mingw_$(WIN)//bin//")

# download
# https://eogdata.mines.edu/nighttime_light/annual/v20/
ntlviirsdir = 'V:/00_GLB/010_imageryBaseMaps/VIIRS/VIIRS2_0'
viirs_tmpdir <- 'D:/temp/viirs/'

require(Rcpp)
today <- '20210301'
  
setwd(wkdir)
output_version <- today


adm0_lake_chad_shp = paste0(wkdir,'/local/gis_data/003_boundaries/lake_chad/lake_chad_adm0')
adm0_reg_shp <- rgdal::readOGR(dirname(adm0_lake_chad_shp), basename(adm0_lake_chad_shp))
adm0_lake_chad_sf <- read_sf(dsn=dirname(adm0_lake_chad_shp), layer=basename(adm0_lake_chad_shp))
fishnet_ext <- extent(adm0_reg_shp)

##
## NTL AS SIZE VAR
##

ntl_type_list = c('median_masked','average_masked','lit_masked','average','median')

for (ntl_type in ntl_type_list){
  tempd <- paste0(viirs_tmpdir ,ntl_type)  
  ntl_list_fn <- list.files(tempd, pattern='.gz', full.names = TRUE)
  ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01')
  names(ntl.zones)
  summary(ntl.zones$OBJECTID)
  countries_GID_0 <- unique(ntl.zones$GID_0)
  # add area km2 #
  ntl.zones$sf_area_km2 <- st_area(ntl.zones) / 1000000
  
  setwd(wkdir)
  require(R.utils)
  
  for (ntl_fn_i in ntl_list_fn) {
    print(ntl_fn_i)
    year_i <- substr(basename(ntl_fn_i),12,15)
    if (!file.exists(tools::file_path_sans_ext(ntl_fn_i))) 
      {r <- crop(raster(R.utils::gunzip(ntl_fn_i, remove = FALSE, temporary= TRUE)),fishnet_ext)} else 
      {r <- crop(raster(tools::file_path_sans_ext(ntl_fn_i)),fishnet_ext)}
    print(cellStats(r, 'min'))
    print('adjust negative values to 0')
    r[r<1]==1
    ntl.zones[paste0("ntlv",year_i)] <- exact_extract(r, ntl.zones, 'mean')
    file.remove(tools::file_path_sans_ext(ntl_fn_i))
    remove(r)
  }
  
  out.df <- as.data.frame(st_drop_geometry(ntl.zones))
  out_fn <- paste0(wkdir,'/proc_data/ntl_viirs/ntl_viirs_grt1_',ntl_type,today,'.dta')
  xfun::dir_create(dirname(out_fn))
  save.dta13(out.df,out_fn)

}